1. Introduzione

Obiettivo del Progetto

Il presente progetto si propone di applicare diverse tecniche di clustering al dataset banknote al fine di identificare gruppi omogenei presenti nei dati e confrontare l’efficacia dei vari metodi. L’analisi intende mettere in pratica le metodologie apprese nel modulo di Data Analysis & Statistical Learning, in particolare, verranno trattati i seguenti algoritmi di clustering: Clustering Gerarchico Agglomerativo, Clustering Partizionale (K-means e K-medoids) e Clustering Model-based (Mixture of Gaussian).

TODO: specificare PCA e mixture of regression

Tra diversi dataset consigliati dai docenti, è stato scelto il dataset banknote, per un motivo puramente professionale, essendo tutti i membri del team legati professionalmente alle banconote. TODO: togliere?

Il documento è organizzato nei seguenti capitoli:

  1. Introduzione: Presentazione del progetto, degli obiettivi e del contesto.
  2. Descrizione del Dataset e Analisi Esplorativa: Esplorazione preliminare dei dati, analisi descrittiva e operazioni di preprocessing.
  3. Applicazione delle Tecniche di Clustering: Implementazione dei metodi di clustering (gerarchico, partizionale e model-based).
  4. Validazione e Visualizzazione dei Risultati: Confronto dei risultati mediante metriche di validazione e visualizzazione dei cluster (con l’ausilio della PCA).
  5. Conclusioni e Sviluppi Futuri: Sintesi dei risultati ottenuti e prospettive per ulteriori analisi.

2. Descrizione del Dataset e Analisi Esplorativa

In questa sezione esamineremo il dataset banknote per comprenderne la struttura, le caratteristiche e le relazioni tra le variabili. Il dataset proviene dal pacchetto mclust e contiene misure che possono essere utilizzate per distinguere tra banconote autentiche e false (o per altri scopi diagnostici). Il dataset contiene le seguenti variabili:

Caricamento del Dataset

library(mclust)
data(banknote)

Esplorazione della Struttura e Statistiche Descrittive del dataset

Visualizziamo le prime righe del dataset

head(banknote)
##    Status Length  Left Right Bottom  Top Diagonal
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4

Esaminiamo la struttura del dataset

str(banknote)
## 'data.frame':    200 obs. of  7 variables:
##  $ Status  : Factor w/ 2 levels "counterfeit",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Length  : num  215 215 215 215 215 ...
##  $ Left    : num  131 130 130 130 130 ...
##  $ Right   : num  131 130 130 130 130 ...
##  $ Bottom  : num  9 8.1 8.7 7.5 10.4 9 7.9 7.2 8.2 9.2 ...
##  $ Top     : num  9.7 9.5 9.6 10.4 7.7 10.1 9.6 10.7 11 10 ...
##  $ Diagonal: num  141 142 142 142 142 ...

Il dataset è composto da 200 osservazioni e 7 variabili, tutte di tipo numerico ad eccezione di ‘Status’ che è una variabile binaria che può assumere uno tra i seguenti valori: ‘counterfeit’ e ‘genuine’.

Visualizziamo un riepilogo statistico delle variabili.

summary(banknote)
##          Status        Length           Left           Right      
##  counterfeit:100   Min.   :213.8   Min.   :129.0   Min.   :129.0  
##  genuine    :100   1st Qu.:214.6   1st Qu.:129.9   1st Qu.:129.7  
##                    Median :214.9   Median :130.2   Median :130.0  
##                    Mean   :214.9   Mean   :130.1   Mean   :130.0  
##                    3rd Qu.:215.1   3rd Qu.:130.4   3rd Qu.:130.2  
##                    Max.   :216.3   Max.   :131.0   Max.   :131.1  
##      Bottom            Top           Diagonal    
##  Min.   : 7.200   Min.   : 7.70   Min.   :137.8  
##  1st Qu.: 8.200   1st Qu.:10.10   1st Qu.:139.5  
##  Median : 9.100   Median :10.60   Median :140.4  
##  Mean   : 9.418   Mean   :10.65   Mean   :140.5  
##  3rd Qu.:10.600   3rd Qu.:11.20   3rd Qu.:141.5  
##  Max.   :12.700   Max.   :12.30   Max.   :142.4

Notiamo che la colonna ‘Status’ è perfettamente bilanciata.

Controlliamo se ci sono valori mancanti nelle colonne.

colSums(is.na(banknote))
##   Status   Length     Left    Right   Bottom      Top Diagonal 
##        0        0        0        0        0        0        0

Come possiamo vedere nessuna colonna presenta valori mancanti.

Analisi Esplorativa dei Dati (EDA)

Questa sezione fornisce una panoramica iniziale del dataset banknote, evidenziando la struttura, le principali statistiche descrittive e le relazioni tra le variabili. Queste informazioni sono fondamentali per orientare le successive fasi di clustering e per capire eventuali necessità di preprocessing (ad es. scaling o trasformazioni) prima dell’applicazione dei metodi di clustering.

Analisi delle Correlazioni

Verifichiamo le correlazioni tra le variabili numeriche per capire se esistono relazioni forti che potrebbero influenzare l’analisi di clustering.

cor_matrix <- cor(banknote[, -1])
print(cor_matrix)
##               Length       Left      Right     Bottom         Top   Diagonal
## Length    1.00000000  0.2312926  0.1517628 -0.1898009 -0.06132141  0.1943015
## Left      0.23129257  1.0000000  0.7432628  0.4137810  0.36234960 -0.5032290
## Right     0.15176280  0.7432628  1.0000000  0.4867577  0.40067021 -0.5164755
## Bottom   -0.18980092  0.4137810  0.4867577  1.0000000  0.14185134 -0.6229827
## Top      -0.06132141  0.3623496  0.4006702  0.1418513  1.00000000 -0.5940446
## Diagonal  0.19430146 -0.5032290 -0.5164755 -0.6229827 -0.59404464  1.0000000
library(reshape2)
library(ggplot2)

cor_melted <- melt(cor_matrix)
ggplot(cor_melted, aes(x = Var1, y = Var2, fill = value)) + geom_tile() +
geom_text(aes(label = round(value, 2)), color = "black",
size = 5) + scale_fill_gradient2(low = "blue", mid = "white",
high = "red", midpoint = 0, limits = c(-1, 1)) + theme_minimal() +
labs(title = "Correlazione variabili", fill = "Correlazione") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))

TODO: Commento:

Osserva i valori di correlazione per identificare variabili fortemente correlate (positivamente o negativamente).

Discuti l’impatto di eventuali correlazioni elevate: potrebbero ridondare in informazioni e influenzare il clustering.

Se trovi correlazioni elevate, suggerisci che una riduzione della dimensionalità (es. tramite PCA) potrebbe aiutare.

Distribuzione delle Variabili

Infine, esaminiamo la distribuzione di ciascuna variabile con degli istogrammi:

# Impostiamo un layout per visualizzare gli istogrammi
par(mfrow = c(2, 2))
for(i in 2:7){
    hist(banknote[[i]], main = paste("Istogramma di", names(banknote)[i]),
    xlab = names(banknote)[i], col = "lightblue", border = "white")
}

par(mfrow = c(1, 1))  # Ripristiniamo il layout originale

TODO: Commento:

Descrivi la forma delle distribuzioni: simmetriche, asimmetriche (skewed), unimodali o multimodali.

Nota eventuali outlier o code lunghe che potrebbero influenzare la successiva analisi di clustering.

Confronta le distribuzioni tra le variabili per capire se alcune variabili hanno una variabilità molto diversa, che potrebbe richiedere una normalizzazione/scaling.

Matrice di Scatterplot

Per avere una prima idea delle relazioni tra le variabili, creiamo una matrice di scatterplot. Utilizziamo la colonna ‘Status’ per colorare i punti.

library(GGally)
ggpairs(banknote[,-1],upper = list(continuous = "density", combo = "box_no_facet"),
        lower = list(continuous = "points", combo = "dot_no_facet"),aes(colour=banknote$Status))

Lungo la diagonale sono visibili le distribuzioni univariate delle variabili. Si nota una sovrapposizione parziale tra i due gruppi di colore (Status), ma anche zone in cui uno dei due colori prevale. Nei grafici di densità bivariata (metà diagonale superiore della matrice) e negli scatter plot (metà diagonale inferiore della matrice), si osservano alcune regioni in cui i punti di uno stesso colore tendono a concentrarsi, suggerendo cluster parzialmente distinti. Nonostante la sovrapposizione in diverse aree, si intravedono regioni in cui un gruppo prevale. Ciò suggerisce che un algoritmo di classificazione o di clustering potrebbe distinguere parzialmente le due classi, anche se non in modo perfetto con una singola coppia di variabili. È importante sottolineare che questo grafico mostra solo le relazioni a coppie. Ciò significa che l’analisi si limita a considerare due variabili alla volta, mentre una clusterizzazione netta potrebbe emergere solo quando si considerano contemporaneamente più variabili. In altre parole, anche se alcune coppie non evidenziano una separazione chiara, combinando più dimensioni si potrebbe scoprire una struttura clusterizzata più marcata. Il fatto che in alcune proiezioni (coppie di variabili) si vedano zone di separazione o contorni distinti suggerisce che esistono sottostrutture nei dati. Notiamo infine che le variabili ‘Bottom’ e ‘Diagonal’ sembrerebbero separare in modo migliore il dataset rispetto alle altre variabili.

D’ora in avanti, per ciascun algoritmo di clustering che proveremo, verranno effettuate tre diverse analisi:

  1. Utilizzando tutte le variabili disponibili nel dataset.

  2. Considerando esclusivamente le variabili ‘Diagonal’ e ‘Bottom’.

  3. Applicando l’Analisi delle Componenti Principali (PCA) per ridurre la dimensionalità.

Queste tre approcci permetteranno di valutare l’impatto delle diverse selezioni di variabili e della riduzione dimensionale sulle prestazioni e sui risultati degli algoritmi di clustering.

Prima di applicare gli algoritmi di clustering, è fondamentale verificare se il dataset presenta una naturale propensione alla formazione di gruppi.

Valutazione della tendenza al clustering del dataset banknote

Utilizziamo specifiche misure statistiche per determinare se i dati sono strutturati in cluster o se, al contrario, sono distribuiti in modo casuale. Questi metodi permettono di valutare se il clustering è effettivamente significativo e utile per l’analisi del dataset.

Hopkins statistics

L’Hopkins Statistic è una misura utilizzata per valutare la tendenza di un dataset alla clusterizzazione, confrontando la sua struttura con quella di un insieme di punti generati casualmente. Assume valori compresi tra 0 e 1. Valori vicini a 0 indicano che il dataset è altamente strutturato in cluster, suggerendo che l’applicazione di algoritmi di clustering può essere efficace. Valori vicini a 0.5 o superiori suggeriscono che i dati sono distribuiti in modo casuale o uniforme, e il clustering potrebbe non essere significativo.

library(clustertend)

banknote_num <- banknote[, -1]

banknote_num <- scale(banknote_num)

set.seed(987)

hopkins_stat <- hopkins(banknote_num, n = nrow(banknote_num) - 1)

random_df <- apply(banknote_num, 2, function(x){runif(length(x), min(x), max(x))})
random_df <- as.data.frame(random_df)

random_df <- scale(random_df)

set.seed(987)

hopkins_stat_rand <- hopkins(random_df, n = nrow(random_df)-1)

cat("Hopkins Statistics on Banknote dataset:", hopkins_stat$H, "\n")
## Hopkins Statistics on Banknote dataset: 0.2642237
cat("Hopkins Statistics on Random dataset:", hopkins_stat_rand$H, "\n")
## Hopkins Statistics on Random dataset: 0.4926061

Un valore di Hopkins pari a 0.2519264, essendo inferiore a 0.5, indica una tendenza alla clusterizzazione. In altre parole, le distanze tra i punti reali sono molto più ridotte rispetto a quelle dei punti generati casualmente, suggerendo la presenza di gruppi ben definiti nel dataset. Questo risultato supporta l’ipotesi che il dataset sia strutturato in cluster distinti.

Visual Assessment of cluster Tendency (VAT)

L’algoritmo VAT produce una matrice di dissimilarità ordinata, ottenuta riorganizzando le righe e le colonne della matrice delle distanze per evidenziare eventuali blocchi di punti “ravvicinati” (ossia potenziali cluster). Visualizziamo la matrice con l’intento di verificare se si formano lungo la diagonale secondaria del grafico, dei blocchi rossi che rappresentano regioni di osservazioni che condividono una forte similarità.

library("factoextra")
library(gridExtra)

p1 <- fviz_dist(dist(banknote), show_labels = FALSE) + labs(title = "Banknote data")
p2 <- fviz_dist(dist(random_df), show_labels = FALSE) + labs(title = "Random data")

grid.arrange(p1, p2, ncol = 2)

Nei grafici VAT il rosso indica alta similarità (ossia bassa dissimilarità), mentre il blu indica bassa similarità (ossia alta dissimilarità). Relativamente al grafico del dataset banknote, notiamo che lungo la diagonale sono presenti regioni rosse che indicano dati simili tra loro, suggerendo la presenza di cluster. La presenza di zone blu che evidenziano una maggiore dissimilarità tra gruppi di dati, supporta l’ipotesi che esistano cluster distinti nel dataset. Inoltre la presenza di transizioni marcate dal rosso al blu conferma che esistono confini ben definiti tra i cluster. Notiamo invece come nel dataset randomico non c’è alcuna struttura a cluster evidente. I dati appaiono disordinati e non presenteranno raggruppamenti chiari, il che indica che un clustering non ha senso o che i dati non sono adatti per essere separati in gruppi distinti.

Possiamo concludere che il dataset Banknote ha una tendezza alla clusterizzazione, procediamo perciò con il calcolo della PCA e successivamente alla clusterizzazione.

PCA

In questa sezione applichiamo la PCA sul dataset, considerando solo le colonne numeriche. Anche se le variabili sono tutte misurate in millimetri e quindi nella stessa scala, applichiamo comunque la standardizzazione per evitare che differenze nella dispersione (variabili con varianze maggiori) influenzino in modo sproporzionato i risultati della PCA.

library("FactoMineR")

banknote_num <- banknote[, -1]
head(banknote_num)
##   Length  Left Right Bottom  Top Diagonal
## 1  214.8 131.0 131.1    9.0  9.7    141.0
## 2  214.6 129.7 129.7    8.1  9.5    141.7
## 3  214.8 129.7 129.7    8.7  9.6    142.2
## 4  214.8 129.7 129.6    7.5 10.4    142.0
## 5  215.0 129.6 129.7   10.4  7.7    141.8
## 6  215.7 130.8 130.5    9.0 10.1    141.4
pca_result <- PCA(banknote_num, scale = TRUE, graph = FALSE)

Visualizzazione e interpretazione dei risultati

Estraiamo gli autovalori e calcoliamo sia la varianza spiegata per ciascuna componente sia la varianza cumulativa. Queste informazioni ci serviranno per decidere il numero ottimale di componenti principali da mantenere.

eigenvalues <- get_eigenvalue(pca_result)[,1]
explained_variance <- eigenvalues / sum(eigenvalues)
cumulative_variance <- cumsum(explained_variance)

print(eigenvalues)
##     Dim.1     Dim.2     Dim.3     Dim.4     Dim.5     Dim.6 
## 2.9455582 1.2780838 0.8690326 0.4497687 0.2686769 0.1888799
print(explained_variance)
##      Dim.1      Dim.2      Dim.3      Dim.4      Dim.5      Dim.6 
## 0.49092637 0.21301396 0.14483876 0.07496145 0.04477948 0.03147998
print(cumulative_variance)
##     Dim.1     Dim.2     Dim.3     Dim.4     Dim.5     Dim.6 
## 0.4909264 0.7039403 0.8487791 0.9237405 0.9685200 1.0000000

Abbiamo ora a disposizione gli autovalori e le informazioni sulla varianza, utili per valutare l’importanza di ciascuna componente principale.

Utilizziamo la Kaiser Rule, una tecnica utilizzata per determinare il numero di componenti principali da mantenere quando si esegue l’analisi delle componenti principali. La regola si basa sull’idea di conservare solo le componenti con una varianza spiegata maggiore di 1. Questo criterio si applica al valore proprio (eigenvalue) delle componenti principali.

kaiser_components <- sum(eigenvalues > 1)
cat("Numero di componenti selezionati con la Kaiser rule:", kaiser_components, "\n")
## Numero di componenti selezionati con la Kaiser rule: 2

Otteniamo 2 come numero di componenti principali da mantenere, ciò significa che le due componenti principali con valori propri superiori a 1 spiegano una quantità di varianza maggiore di quella che verrebbe spiegata da una singola variabile originale del dataset.

Un ulteriore criterio è quello di selezionare il numero minimo di componenti la cui somma di varianza spiegata (varianza comulativa), raggiunge almeno l’80% del totale. In questo modo garantiamo una buona rappresentazione dei dati.

variance_threshold <- min(which(cumulative_variance >= 0.80))

cat("Numero di componenti principali che spiegano l'80% della varianza comulativa:", variance_threshold, "\n")
## Numero di componenti principali che spiegano l'80% della varianza comulativa: 3
fviz_eig(pca_result, addlabels = TRUE, ylim = c(0, 50))

Nel grafico vediamo la varianza spiegata da ciascuna componente. Con le prime tre componenti principali riusciamo a spiegare l’84.9% della varianza cumulativa.

Vediamo infine un ultimo metodo per scegliere il numero di componenti principali: l’elbow method. Mediante un grafico scree plot visualizziamo gli autovalori per ciascuna componente e individuiamo “il gomito”, ovvero il punto nel grafico in cui la diminuzione degli autovalori rallenta. Tale punto può essere interpretato come un’indicazione del numero ottimale di componenti. Includiamo nel grafico anche due linee verticali per evidenziare i criteri della varianza cumulativa (rosso) e della Kaiser Rule (blu).

scree_plot <- data.frame(PC = 1:length(eigenvalues), Eigenvalue = eigenvalues)

ggplot(scree_plot, aes(x = PC, y = Eigenvalue)) +
    geom_point() +
    geom_line() +
    geom_vline(xintercept = variance_threshold, linetype = "dashed", color = "red") +
    geom_vline(xintercept = kaiser_components, linetype = "dotted", color = "blue") +
    labs(title = "Scree Plot for PCA", x = "Principal Component", y = "Eigenvalue") +
    theme_minimal()

Nel grafico il gomito sembra collocarsi intorno alla quarta componente. I tre criteri danno perciò risultati differenti, per mantenere un equilibrio tra interpretabilità e preservazione della varianza, si potrebbe optare per 3 componenti.
Di seguito viene costruito il nuovo dataset in cui le variabili sono sostituite con le prime 3 PCA, che verrà usato più avanti nelle analisi.

pca_data <- as.data.frame(pca_result$ind$coord[, 1:3])
colnames(pca_data) <- c("PCA1", "PCA2", "PCA3")
banknote_pca <- cbind(pca_data, Status = banknote$Status)
banknote_pca <- banknote_pca %>% dplyr::select(Status, everything())
head(banknote_pca)
##    Status       PCA1        PCA2       PCA3
## 1 genuine  1.7474012  1.65082829  1.4237611
## 2 genuine -2.2743177 -0.53879328  0.5326484
## 3 genuine -2.2774015 -0.10767707  0.7174149
## 4 genuine -2.2835546 -0.08765431 -0.6056336
## 5 genuine -2.6321283  0.03919590  3.1963847
## 6 genuine  0.7584073  3.08874513  0.7864804

Per comprendere meglio il contributo delle variabili alle componenti principali, visualizziamo le loro coordinate e realizziamo una mappa di correlazione. Questo ci aiuta a capire come le variabili si relazionano alle componenti principali.

var <- get_pca_var(pca_result)

var$coord
##                Dim.1      Dim.2       Dim.3       Dim.4      Dim.5
## Length   -0.01199158  0.9219364 -0.01648225 -0.38536590  0.0304764
## Left      0.80279596  0.3866019  0.09637548  0.26485400 -0.3314768
## Right     0.83526859  0.2854104  0.11510550  0.28856524  0.3183114
## Bottom    0.69810421 -0.3009779  0.54398559 -0.27072283  0.1116898
## Top       0.63139786 -0.1034278 -0.73418921 -0.07392333  0.1139569
## Diagonal -0.84690418  0.3096965  0.10615680  0.26284740  0.1763188
fviz_pca_var(pca_result, col.var = "contrib",
            gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

Questo passaggio ci offre una visione immediata di come ogni variabile contribuisce alle componenti principali e della loro interrelazione.

Approfondiamo l’analisi mostrando i contributi delle variabili attraverso diversi metodi grafici, in modo da identificare quali variabili influenzano maggiormente ciascuna componente.

print(round(var$contrib,digits = 3))
##           Dim.1  Dim.2  Dim.3  Dim.4  Dim.5
## Length    0.005 66.503  0.031 33.019  0.346
## Left     21.880 11.694  1.069 15.596 40.896
## Right    23.686  6.374  1.525 18.514 37.712
## Bottom   16.545  7.088 34.052 16.295  4.643
## Top      13.534  0.837 62.027  1.215  4.833
## Diagonal 24.350  7.504  1.297 15.361 11.571
library("corrplot")
corrplot(var$contrib, is.corr=FALSE)

L’analisi dei contributi evidenzia il peso relativo di ciascuna variabile nelle componenti principali, utile per interpretare i risultati della PCA.

Visualizziamo i contributi delle variabili specificamente per la prima componente principale.

fviz_contrib(pca_result, choice = "var", axes = 1)

Similmente al passaggio precedente, visualizziamo i contributi delle variabili per la seconda componente principale.

fviz_contrib(pca_result, choice = "var", axes = 2)

3. Applicazione delle Tecniche di Clustering

3.1 Clustering Gerarchico Agglomerativo

Il Clustering Gerarchico Agglomerativo (HAC) è una tecnica di clustering che costruisce una gerarchia di cluster attraverso un approccio “bottom-up”. Inizia considerando ogni punto dati come un cluster separato e, successivamente, unisce iterativamente i cluster più simili fino a formare un unico cluster che racchiude tutti i dati. Questo metodo non richiede la specifica del numero di cluster a priori e produce una rappresentazione ad albero, nota come dendrogramma, che facilita la visualizzazione delle relazioni tra i cluster.

Per applicare il clustering gerarchico agglomerativo al nostro dataset, seguiamo questi passaggi:

  1. Scaling:

  2. Calcolo della Matrice di Distanza: Determiniamo le distanze tra tutte le coppie di punti nel dataset. Una misura comune è la distanza euclidea.

  3. Applicazione del clustering agglomerativo: Utilizziamo la funzione hclust() per eseguire il clustering gerarchico e generare il dendrogramma.

  4. Taglio del Dendrogramma: Decidiamo il numero ottimale di cluster utilizzando varie tecniche, e tagliamo il dendrogramma al livello appropriato.

  5. Valutazione: Valutiamo i risultati del clustering utilizzando varie metriche.

Clustering con tutte le variabili

Scaliamo il dataset e calcoliamo la relativa matrice di dissimilatità mediante distanza euclidea e la mostriamo. si è scelta la distanza euclidea poichè tutte le variabili del dataset su cui applicheremo il clustering sono di tipo numerico.

library(pheatmap)

df <- scale(banknote_num)

distanze <- dist(df, method = "euclidean")

diss_matrix <- as.matrix(distanze)

pheatmap(diss_matrix,
         cluster_rows = FALSE,
         cluster_cols = FALSE,
         show_rownames = FALSE,
         show_colnames = FALSE,
         color = colorRampPalette(c("blue", "white", "red"))(50),
         border_color = NA)

Applichiamo il clustering con metodo di linkage ‘ward’ che generalmente insieme all’average linkage, performa meglio.

hclust_clustering <- hclust(distanze, method = "ward.D2")

Guardiamo il dendogramma cercando di identificare il numero ottimale di cluster. Utilizzeremo poi anche altri metodi che ci premetteranno di prendere una scelta migliore.

Numero di cluster

fviz_dend(hclust_clustering, show_labels = FALSE, as.ggplot = TRUE)

Il dendrogramma risultante mostra come le osservazioni vengono unite progressivamente in cluster. L’altezza dei rami indica la distanza a cui avviene la fusione: rami più bassi corrispondono a fusioni tra punti o cluster molto simili. Analizzando il dendrogramma, possiamo decidere il numero di cluster tagliando l’albero a un’altezza che separa i gruppi in modo significativo. Dall’osservazione del dendrogramma, si nota un salto particolarmente ampio a un’altezza intorno a 30, il che indica che fino a quella soglia i gruppi sono relativamente compatti e ben separati. Tagliare il dendrogramma a quel livello produrrebbe 2 grandi cluster. Volendo una suddivisione più fine, si può individuare un altro salto meno marcato (intorno a 15) che porterebbe a 3 cluster. Procediamo con l’analisi del numero ottimale di cluster utilizzando altri metodi.

Osserviamo l’andamento del Total Within Sum of Squares (WSS) rispetto a k, cercando di individuare il punto in cui la curva inizia a “piegarsi” più marcatamente, cioè dove l’aggiunta di ulteriori cluster non porta a una riduzione significativa del WSS. Di seguito il grafico Elbow.

fviz_nbclust(df, FUN = hcut, method = "wss") +
    geom_vline(xintercept = 3, linetype = 2) +
    labs(subtitle = "Elbow method")

Tra k = 1 e k = 2, la diminuzione del WSS è molto forte (ciò è normale, perché si passa da un unico cluster a due). Tra k = 2 e k = 3, c’è ancora un calo netto del WSS. Oltre k = 3, la pendenza inizia a ridursi gradualmente. Dal grafico, sembra che il calo principale avvenga per k = 3, mentre oltre si riduce più lentamente. Per k = 3 si ha un notevole miglioramento rispetto a 2, ma passando a 4 (e oltre) il guadagno in termini di riduzione del WSS diventa più contenuto.

Analizziamo ora il grafico dell’Average Silhouette Width.

fviz_nbclust(df, FUN = hcut, method = "silhouette") +
    labs(subtitle = "Silhouette method")

L’indice di silhouette risulta più alto in corrispondenza di k = 2. Questo suggerisce che la separazione tra i cluster e la coesione interna sia massima quando si dividono i dati in due gruppi. Con k = 3, l’indice diminuisce, e oltre tale valore rimane su livelli inferiori (intorno a 0.2-0.3), senza tornare ai valori di picco.

Osserviamo il grafico del Gap Statistic che confronta la dispersione all’interno dei cluster ottenuta dal clustering, con quella attesa in un dataset generato casualmente (quindi privo di gruppi). Ossia calcola la differenza (gap) tra il log della dispersione osservata e quella media attesa sotto una distribuzione nulla.

set.seed(123)
fviz_nbclust(df, FUN = hcut, method = "gap_stat", nboot = 500)+
    labs(subtitle = "Gap statistic method")

Nel grafico, la linea che rappresenta il punto in cui il Gap Statistic non aumenta più in maniera significativa, è su 3 cluster, il che indica che quando i dati vengono divisi in 3 gruppi la struttura di clustering risulta migliore rispetto a quella attesa in un dataset casuale. In altre parole, il modello a 3 cluster spiega meglio la struttura intrinseca dei dati.

Utilizziamo infine il pacchetto ‘NbClust’ che valuta il numero ottimale di cluster utilizzando numerosi indici interni. Tra questi, alcuni indici come l’indice di Hubert e il Dindex vengono visualizzati graficamente.

library("NbClust")

nb <- NbClust(df, distance = "euclidean", min.nc = 2, max.nc= 10, method = "ward.D2")

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 7 proposed 2 as the best number of clusters 
## * 10 proposed 3 as the best number of clusters 
## * 1 proposed 4 as the best number of clusters 
## * 1 proposed 7 as the best number of clusters 
## * 1 proposed 8 as the best number of clusters 
## * 3 proposed 10 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  3 
##  
##  
## *******************************************************************

La decisione finale basata sulla majority rule suggerisce 3 cluster, indicando che, complessivamente, questa soluzione offre una migliore separazione e coesione interna rispetto alle altre soluzioni.

Complessivamente la maggioranza dei metodi analizzati (elbow, gap statistic, nbclust) converge su una soluzione a 3 cluster. Pertanto optiamo per una scelta di “maggioranza” e scegliamo perciò un numero di cluster pari a 3.

Taglio del dendrogramma

Procediamo con il taglio del dendogramma in modo da ottenere 3 cluster. Successivamente visualizziamo una tabella di frequenza che mostra il numero di osservazioni assegnate a ciascun cluster, permettendoci di verificare la distribuzione dei dati tra i 3 cluster ottenuti.

cluster_cut <- cutree(hclust_clustering, k = 3)

table(cluster_cut)
## cluster_cut
##   1   2   3 
## 102  76  22
fviz_dend(hclust_clustering, k = 3,
        cex = 0.5,
        k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
        color_labels_by_k = TRUE,
        rect = TRUE
)

La maggior parte dei dati è stata assegnata ai primi due cluster, mentre il terzo cluster contiene un numero significativamente minore di osservazioni. Questo potrebbe suggerire che il cluster 3 rappresenta un gruppo più piccolo, con caratteristiche potenzialmente diverse o più specifiche rispetto ai gruppi più grandi.

Vediamo i risultati di clustering nello spazio originale, tramite la matrice di scatterplot a coppie con i punti colorati in base ai cluster di appartenenza.

ggpairs(df, aes(color = as.factor(cluster_cut))) +
  ggtitle("Matrice Pairwise Scatterplot con Cluster")

Valutazione

Valutiamo ora la qualità dei risultati del nostro algoritmo di clustering mediante quattro metriche: Average Silhouette Width, Dunn Index, Corrected Rand Index e Meila’s Variation of Information Index.

Iniziamo con la metrica della silhouette media (calcolata come la media di s(i) su tutte le osservazioni), che fornisce un’indicazione complessiva della qualità del clustering: si auspica un valore il più vicino ad 1 possibile.

library(cluster)

sil <- silhouette(cluster_cut, distanze)

fviz_silhouette(sil, palette = "jco", ggtheme = theme_classic())
##   cluster size ave.sil.width
## 1       1  102          0.30
## 2       2   76          0.30
## 3       3   22          0.39

Un valore di 0.31 suggerisce che la struttura dei cluster non è particolarmente forte: c’è una separazione moderata, ma c’è anche una certa sovrapposizione tra i cluster.

Analizziamo ora il Dunn Index, che valuta la separazione tra i cluster rispetto alla loro compattezza, cioè quanto i cluster siano ben separati e internamente compatti. Se i dati contengono cluster compatti e ben separati, ci si aspetta che il diametro dei cluster sia piccolo e che la distanza tra i cluster sia grande. Pertanto l’indice di Dunn deve essere massimizzato.

library(clusterCrit)

dunn_index <- intCriteria(as.matrix(df), cluster_cut, "Dunn")$dunn

print(paste("Dunn Index:", dunn_index))
## [1] "Dunn Index: 0.1250576993898"

Un Dunn Index di 0.125 indica che i cluster hanno una separazione piuttosto bassa rispetto alla loro coesione interna. Una bassa separazione tra cluster suggerisce possibili problemi di sovrapposizione.

Analizziamo il Corrected Rand Index che misura quanto bene il clustering ottenuto corrisponde a una partizione di riferimento. Nel nostro caso utilizziamo la colonna ‘Status’ del dataset. Ed infine il Meila’s Variation of Information (VI) Index che misura quanto due cluster differiscono tra loro in termini di entropia. A differenza del Corrected Rand Index, che valuta la similarità, il VI Index misura la distanza tra due partizioni (va perciò minimizzato).

library("fpc")

status <- as.numeric(banknote$Status)

clust_stats <- cluster.stats(d = dist(df), status, cluster_cut)

print(paste("Corrected Rand Index:", clust_stats$corrected.rand))
## [1] "Corrected Rand Index: 0.791980314729782"
print(paste("Meila's Index:", clust_stats$vi))
## [1] "Meila's Index: 0.359179851985559"

Un Corrected Rand Index di 0.792 indica che il clustering è molto simile alla vera classificazione dei dati. Un Meila’s Variation of Information Index di 0.359 è un valore relativamente basso, il che significa che i cluster trovati non sono troppo distanti dalla struttura effettiva dei dati. Il fatto che il Corrected Rand Index sia alto (0.792) e il VI Index sia basso conferma che il clustering ha una buona coerenza con la partizione reale. Tuttavia, poiché Dunn Index e Silhouette suggeriscono che la separazione tra cluster è debole, questo potrebbe indicare sovrapposizione tra i cluster, ma senza compromettere eccessivamente la qualità dell’assegnazione.

Ottimizzazione

Ottimizziamo il clustering esplorando diverse combinazioni di metodi di linkage e numeri di cluster. Per ciascuna combinazione, calcoliamo le quattro metriche di valutazione precedentemente analizzate, ottenendo un dataframe con tutti i risultati. Successivamente, identifichiamo il metodo di clustering migliore per ciascuna metrica. Poiché le diverse metriche potrebbero suggerire metodi differenti, assegniamo un ranking a ciascuna di esse e calcoliamo un ranking totale sommando i punteggi (invertendo il criterio per le metriche da massimizzare). Infine, ordiniamo i risultati e selezioniamo la combinazione con il punteggio complessivo più basso, individuando così il miglior clustering in base a tutte le metriche disponibili.

library(dplyr)

evaluate_clustering <- function(df, distanze, status, linkage_methods = c("ward.D2", "single", "complete", "average", "centroid"), min_k = 2, max_k = 10) {
    results_all <- data.frame(Method = character(),
                            k = numeric(),
                            Silhouette = numeric(),
                            Dunn = numeric(),
                            Corrected_Rand = numeric(),
                            VI = numeric(),
                            stringsAsFactors = FALSE)

    for(method in linkage_methods) {
        hclust_obj <- hclust(distanze, method = method)
        
        for(k in min_k:max_k) {
            cluster_cut <- cutree(hclust_obj, k = k)
            
            sil <- silhouette(cluster_cut, distanze)
            avg_sil <- mean(sil[, 3])
            
            dunn_index <- intCriteria(as.matrix(df), cluster_cut, "Dunn")$dunn
            
            clust_stats <- cluster.stats(distanze, status, cluster_cut)

            corrected_rand <- clust_stats$corrected.rand

            vi_val <- clust_stats$vi
            
            results_all <- rbind(results_all, data.frame(Method = method,
                                                        k = k,
                                                        Silhouette = avg_sil,
                                                        Dunn = dunn_index,
                                                        Corrected_Rand = corrected_rand,
                                                        VI = vi_val,
                                                        stringsAsFactors = FALSE))
        }
    }
    return(results_all)
}

results_all <- evaluate_clustering(df, distanze, status)

best_silhouette <- results_all %>% group_by(Method) %>% filter(Silhouette == max(Silhouette))
print(best_silhouette)
## # A tibble: 5 × 6
## # Groups:   Method [5]
##   Method       k Silhouette  Dunn Corrected_Rand     VI
##   <chr>    <int>      <dbl> <dbl>          <dbl>  <dbl>
## 1 ward.D2      2      0.375 0.211          0.960 0.0982
## 2 single       2      0.333 0.341          0     0.718 
## 3 complete     2      0.346 0.125          0.591 0.566 
## 4 average      2      0.333 0.341          0     0.718 
## 5 centroid     2      0.333 0.341          0     0.718
best_dunn <- results_all %>% group_by(Method) %>% filter(Dunn == max(Dunn))
print(best_dunn)
## # A tibble: 5 × 6
## # Groups:   Method [5]
##   Method       k Silhouette  Dunn Corrected_Rand     VI
##   <chr>    <int>      <dbl> <dbl>          <dbl>  <dbl>
## 1 ward.D2      2      0.375 0.211          0.960 0.0982
## 2 single       2      0.333 0.341          0     0.718 
## 3 complete    10      0.186 0.179          0.268 1.45  
## 4 average      2      0.333 0.341          0     0.718 
## 5 centroid     2      0.333 0.341          0     0.718
best_corrected_rand <- results_all %>% group_by(Method) %>% filter(Corrected_Rand == max(Corrected_Rand))
print(best_corrected_rand)
## # A tibble: 5 × 6
## # Groups:   Method [5]
##   Method       k Silhouette  Dunn Corrected_Rand     VI
##   <chr>    <int>      <dbl> <dbl>          <dbl>  <dbl>
## 1 ward.D2      2     0.375  0.211       0.960    0.0982
## 2 single       8    -0.0124 0.225       0.000320 0.896 
## 3 complete     3     0.327  0.138       0.688    0.596 
## 4 average      5     0.316  0.222       0.951    0.140 
## 5 centroid     3     0.228  0.298       0.000101 0.742
best_vi <- results_all %>% group_by(Method) %>% filter(VI == min(VI))
print(best_vi)
## # A tibble: 5 × 6
## # Groups:   Method [5]
##   Method       k Silhouette  Dunn Corrected_Rand     VI
##   <chr>    <int>      <dbl> <dbl>          <dbl>  <dbl>
## 1 ward.D2      2      0.375 0.211          0.960 0.0982
## 2 single       2      0.333 0.341          0     0.718 
## 3 complete     2      0.346 0.125          0.591 0.566 
## 4 average      5      0.316 0.222          0.951 0.140 
## 5 centroid     2      0.333 0.341          0     0.718
library(dplyr)

results_ranked <- results_all %>%
  mutate(
    rank_sil = rank(-Silhouette, ties.method = "min"),
    rank_dunn = rank(-Dunn, ties.method = "min"),
    rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
    
    rank_vi = rank(VI, ties.method = "min"),
    
    total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
  )

best_overall <- results_ranked %>% arrange(total_rank) %>% head(1)
print(best_overall)
##    Method k Silhouette      Dunn Corrected_Rand         VI rank_sil rank_dunn
## 1 ward.D2 2  0.3748428 0.2107915      0.9602001 0.09823913        1        21
##   rank_corrected_rand rank_vi total_rank
## 1                   1       1         24

Il risultato mostra che il clustering ottimale secondo queste metriche è Ward.D2 con 2 cluster, poiché bilancia al meglio la coesione, la separazione e la corrispondenza con la verità di riferimento. La combinazione scelta restituisce i valori migliori possibili per le metriche Silhouette, Corrected Rand Index e VI (Rank 1), mentre risulta meno performante rispetto ad altri metodi, ma compensato dagli altri punteggi, per il Dunn Index (Rank 21).

Effettuiamo il clustering con la combinazione risultante e visualizziamo mediante il grafico i cluster ottenuti.

final <- eclust(df, "hclust", k = 2, hc_metric = "euclidean", hc_method = "ward.D2", graph = FALSE)

fviz_cluster(final, geom = "point", ellipse.type = "norm", palette="jco", ggtheme = theme_minimal())

Clustering con PCA

Vediamo come cambia il clustering se invece che considerare tutte le variabili del dataset, consideriamo solo le tre componenti principali trovate prima.

banknote_pca_num <- banknote_pca[,-1]

results_all_pca <- evaluate_clustering(banknote_pca_num, dist(banknote_pca_num), status)

results_ranked_pca <- results_all_pca %>%
  mutate(
    rank_sil = rank(-Silhouette, ties.method = "min"),
    rank_dunn = rank(-Dunn, ties.method = "min"),
    rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
    rank_vi = rank(VI, ties.method = "min"),
    
    total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
  )

best_overall_pca <- results_ranked_pca %>% arrange(total_rank) %>% head(1)
print(best_overall_pca)
##    Method k Silhouette      Dunn Corrected_Rand        VI rank_sil rank_dunn
## 1 average 2   0.431344 0.1056305      0.9020088 0.1997658        1        27
##   rank_corrected_rand rank_vi total_rank
## 1                   2       1         31

Il risultato indica che, applicando il clustering agglomerativo con numero di cluster pari a 2 e metodo di linkage “average” al dataset ridotto a tre componenti principali, otteniamo la combinazione migliore in base ai 4 indici di valutazione, anche se i risultati presentano alcune peculiarità. Rispetto alle altre combinazioni testate, la silhouette e la VI hanno ottenuto i migliori punteggi (rank 1), il Corrected Rand è quasi il migliore (rank 2), invece l’indice Dunn, in questa combinazione, non performa bene rispetto ad altre configurazioni testate (rank 27). Rispetto al clustering agglomerativo applicato su tutte le variabili del dataset il clustering sul dataset ridotto mostra: - Un indice silhouette leggermente migliore (0.431) rispetto a quello su tutte le variabili (0.375), suggerendo che i punti sono, mediamente, meglio raggruppati all’interno dei cluster nel caso PCA. - Il clustering su tutte le variabili ha un indice Dunn migliore (0.211 vs 0.106), il che significa che la separazione minima tra cluster è maggiore e i cluster risultano più distinti. - Entrambe le soluzioni mostrano valori molto alti del Corrected Rand e bassi del VI, ma la soluzione con tutte le variabili è leggermente superiore (0.960 vs 0.902 e 0.098 vs 0.200), indicando una migliore aderenza ad una struttura di riferimento o una maggiore stabilità del clustering.

Mentre l’approccio con PCA (riduzione a 3 componenti) offre cluster con una buona compattezza (indice silhouette migliore), il clustering basato su tutte le variabili, utilizzando il metodo ward.D2, risulta complessivamente superiore per quanto riguarda la separazione tra cluster e la corrispondenza con una struttura di riferimento (migliori valori di Dunn, Corrected Rand e VI).

Vediamo graficamente come sono separati i due cluster.

library(scatterplot3d)

final_pca <- eclust(banknote_pca_num, "hclust", k = 2, hc_metric = "euclidean", hc_method = "average", graph = FALSE)

clusters <- final_pca$cluster

scatterplot3d(banknote_pca_num,
              color = clusters,
              pch = 19,
              main = "Clustering su PCA (PC1, PC2, PC3)",
              xlab = "PC1", ylab = "PC2", zlab = "PC3")

legend("topright", legend = unique(clusters), col = unique(clusters), pch = 19, title = "Cluster")

Clustering con solo le variabili Bottom e Diagonal

Vediamo infine come cambia il clustering considerando solo le variabili ‘Bottom’ e ‘Diagonal’ che sembravano suddividere meglio i dati in gruppi.

banknote_bottom_diagonal <- banknote %>% select(Bottom, Diagonal)

banknote_bd_scaled <- scale(banknote_bottom_diagonal)

results_all_bd <- evaluate_clustering(banknote_bd_scaled, dist(banknote_bd_scaled), status)

results_ranked_bd <- results_all_bd %>%
  mutate(
    rank_sil = rank(-Silhouette, ties.method = "min"),
    rank_dunn = rank(-Dunn, ties.method = "min"),
    rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
    rank_vi = rank(VI, ties.method = "min"),
    
    total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
  )

best_overall_bd <- results_ranked_bd %>% arrange(total_rank) %>% head(1)
print(best_overall_bd)
##    Method k Silhouette       Dunn Corrected_Rand        VI rank_sil rank_dunn
## 1 ward.D2 2  0.6245096 0.09556396         0.9602 0.1120031        1        11
##   rank_corrected_rand rank_vi total_rank
## 1                   2       2         16

In questo caso la combinazione migliore è data da metodo di linkage ‘ward.D2’ e numero cluster pari a 2. Notiamo che in questo caso che il valore Silhouette 0.6245 è più alto rispetto agli altri due approcci (0.431 per PCA e 0.375 usando tutte le variabili), il che indica che i punti all’interno dei cluster sono ben raggruppati e distinti l’uno dall’altro. Ciò supporta l’idea iniziale che queste due variabili siano in grado di dividere bene il dataset. Il Dunn index risulta più basso rispetto a quello ottenuto con tutte le variabili (0.211) o con la PCA (0.106), il che potrebbe significare che alcuni cluster risultano più vicini o con una certa dispersione interna. Il Corrected Rand 0.9602 rimane molto elevato, segnalando una forte coerenza con la struttura di riferimento attesa. L’indice VI è basso, confermando una buona similarità o una discreta aderenza a una struttura di riferimento.

Visualizziamo come il dataset è stato diviso in cluster.

final_bd <- eclust(banknote_bd_scaled, "hclust", k = 2, hc_metric = "euclidean", hc_method = "ward.D2", graph = FALSE)

fviz_cluster(final_bd, geom = "point", ellipse.type = "norm", palette="jco", ggtheme = theme_minimal())

3.2 Clustering Gerarchico Divisivo

Il clustering divisivo è una tecnica gerarchica “top-down” che parte da un unico cluster contenente l’intero dataset e procede a suddividerlo progressivamente in cluster più piccoli. A differenza del clustering agglomerativo, che inizia con ciascun punto come cluster individuale e li unisce gradualmente, il metodo divisivo si concentra sulla separazione iterativa dei dati in gruppi distinti. Ad ogni suddivisione, si cerca di massimizzare la differenza (o la distanza) tra i cluster risultanti e minimizzare la varianza all’interno di ciascun cluster. Il risultato finale è una struttura ad albero (il dendrogramma) che rappresenta le varie divisioni e il livello di dettaglio raggiunto, offrendo una visione della struttura gerarchica dei dati.

Clustering Gerarchico Divisivo con tutte le variabili

Applichiamo il clusering gerarchico divisivo sul nostro dataset Banknote con tutte le variabili. Utilizziamo l’algoritmo DIANA (DIvisive Analysis Clustering) sul dataset standardizzato (scaling) in precedenza e distanza euclidea per calcolare la distanza tra le osservazioni. La funzione ‘eclust’, se non si specifica il parametro k, calcola automaticamente la statistica del gap per stimare il numero ottimale di cluster. Fornisce inoltre informazioni sulla silhouette per tutti i metodi di partizionamento e per il clustering gerarchico.

diana_clustering <- eclust(df, "diana", stand = FALSE, hc_metric = "euclidean")

fviz_dend(diana_clustering, cex = 0.6, show_labels = FALSE, as.ggplot = TRUE)

Osserviamo che la funzione eclust ha stabilito che il numero ottimale di cluster è 3. Proviamo a calcolare la migliore combinazione di numero di cluster con una leggera modifica del codice visto prima.

evaluate_diana_clustering <- function(df, distanze, status, min_k = 2, max_k = 10) {
    results_all <- data.frame(  k = numeric(),
                                Silhouette = numeric(),
                                Dunn = numeric(),
                                Corrected_Rand = numeric(),
                                VI = numeric(),
                                stringsAsFactors = FALSE)
    
    diana_obj <- diana(distanze)

    hclust_obj <- as.hclust(diana_obj)
    
    for(k in min_k:max_k) {
        cluster_cut <- cutree(hclust_obj, k = k)
        
        sil <- silhouette(cluster_cut, distanze)
        avg_sil <- mean(sil[, 3])

        dunn_index <- intCriteria(as.matrix(df), cluster_cut, "Dunn")$dunn
        
        clust_stats <- cluster.stats(distanze, status, cluster_cut)
        corrected_rand <- clust_stats$corrected.rand
        vi_val <- clust_stats$vi
        
        results_all <- rbind(results_all, data.frame(k = k,
                                                     Silhouette = avg_sil,
                                                     Dunn = dunn_index,
                                                     Corrected_Rand = corrected_rand,
                                                     VI = vi_val,
                                                     stringsAsFactors = FALSE))
    }
    return(results_all)
}

results_all_diana <- evaluate_diana_clustering(df, distanze, status)

results_all_diana <- results_all_diana %>%
  mutate(
    rank_sil = rank(-Silhouette, ties.method = "min"),
    rank_dunn = rank(-Dunn, ties.method = "min"),
    rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
    rank_vi = rank(VI, ties.method = "min"),
    
    total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
  )

best_overall_diana <- results_all_diana %>% arrange(total_rank) %>% head(1)
print(best_overall_diana)
##   k Silhouette      Dunn Corrected_Rand        VI rank_sil rank_dunn
## 1 2  0.3723112 0.1819793      0.9406015 0.1540906        1         1
##   rank_corrected_rand rank_vi total_rank
## 1                   1       1          4

Anche in questo caso il risultato suggerisce che il miglior numero di cluster del il dataset è 2. I ranghi di Silhouette, Dunn, Corrected Rand e VI sono tutti pari a 1, il che indica che DIANA con k = 2 ha ottenuto i migliori risultati per queste metriche rispetto agli altri valori di k esplorati. Con un clustering gerarchico divisivo su tutte le variabili numeriche del dataset, con numero di cluster pari a 2 otteniamo: - Silhouette pari a 0.3723. Valore moderato, indica una separazione discreta tra i cluster. - Indice di Dunn pari a 0.18198. Valore piuttosto basso, suggerisce che la separazione tra i cluster non è eccezionale. - Corrected Rand Index pari a 0.9406. Valore alto, indica una buona concordanza con la partizione ideale. - Variational Information (VI) pari a 0.1541. Valore basso, suggerisce una buona corrispondenza con la partizione di riferimento.

Visualizziamo come il dataset è stato diviso in cluster.

final_diana_clustering <- eclust(df, "diana", stand = FALSE, k = 2, hc_metric = "euclidean")

fviz_cluster(final_diana_clustering, geom = "point", ellipse.type = "norm", palette="jco", ggtheme = theme_minimal())

Clustering Gerarchico Divisivo con le componenti principali

Applichiamo il clusering gerarchico divisivo sul nostro dataset considerando solo le 3 componenti principali.

results_pca_diana <- evaluate_diana_clustering(banknote_pca_num, dist(banknote_pca_num), status)

results_pca_diana <- results_pca_diana %>%
  mutate(
    rank_sil = rank(-Silhouette, ties.method = "min"),
    rank_dunn = rank(-Dunn, ties.method = "min"),
    rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
    rank_vi = rank(VI, ties.method = "min"),
    
    total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
  )

best_pca_diana <- results_pca_diana %>% arrange(total_rank) %>% head(1)
print(best_pca_diana)
##   k Silhouette       Dunn Corrected_Rand        VI rank_sil rank_dunn
## 1 2  0.4393924 0.09331653      0.8456292 0.2819728        1         1
##   rank_corrected_rand rank_vi total_rank
## 1                   1       1          4

Ancora una volta viene indicato 2 come numero ideali di cluster.

I risultati ottenuti applicando PCA prima di eseguire il clustering DIANA per k = 2 sono i seguenti:

  • Silhouette pari a 0.4394, leggermente migliore rispetto a quando non è stata applicata la PCA, suggerendo una separazione più chiara tra i cluster.

  • Indice di Dunn pari a 0.0933, piuttosto basso, indicante una separazione limitata tra i cluster.

  • Corrected Rand Index pari a 0.8456, buona concordanza con la partizione ideale, ma inferiore a quando non si applica la PCA.

  • Variational Information (VI) pari a 0.2820, un po’ più alta rispetto al caso senza PCA, indicando una leggera differenza rispetto alla partizione di riferimento.

L’applicazione della PCA prima del clustering DIANA porta a un miglioramento nella separazione dei cluster, come evidenziato dal valore di Silhouette più alto (0.4394 rispetto a 0.3723), ma l’indice di Dunn rimane basso, suggerendo che i cluster non sono completamente separati. Il Corrected Rand Index e il Variational Information (VI) mostrano una discreta corrispondenza con una partizione ideale, ma comunque inferiore rispetto ai risultati ottenuti senza PCA.

Visualizziamo come il dataset è stato diviso in cluster.

final_diana_clustering_pca <- eclust(banknote_pca_num, "diana", stand = FALSE, k = 2, hc_metric = "euclidean")

clusters <- final_diana_clustering_pca$cluster

scatterplot3d(banknote_pca_num,
              color = clusters,
              pch = 19,
              main = "DIANA Clustering su PCA (PC1, PC2, PC3)",
              xlab = "PC1", ylab = "PC2", zlab = "PC3")

legend("topright", legend = unique(clusters), col = unique(clusters), pch = 19, title = "Cluster")

Clustering Gerarchico Divisivo con dataset ridotto

Applichiamo il clusering gerarchico divisivo sul nostro dataset ridotto alle sole due variabili ‘Bottom’ e ‘Diagonal’.

results_bd_diana <- evaluate_diana_clustering(banknote_bd_scaled, dist(banknote_bd_scaled), status)

results_bd_diana <- results_bd_diana %>%
  mutate(
    rank_sil = rank(-Silhouette, ties.method = "min"),
    rank_dunn = rank(-Dunn, ties.method = "min"),
    rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
    rank_vi = rank(VI, ties.method = "min"),
    
    total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
  )

best_bd_diana <- results_bd_diana %>% arrange(total_rank) %>% head(1)
print(best_bd_diana)
##   k Silhouette       Dunn Corrected_Rand        VI rank_sil rank_dunn
## 1 2  0.6241219 0.07897225      0.9406015 0.1540906        2         2
##   rank_corrected_rand rank_vi total_rank
## 1                   1       1          6

I risultati ottenuti applicando DIANA al dataset composto solo dalle variabili ‘Bottom’ e ‘Diagonal’ per k = 2 sono i seguenti:

  • Silhouette pari a 0.6241, valore buono, indica una separazione più marcata tra i cluster rispetto ai precedenti risultati.

  • Indice di Dunn pari a 0.07897, valore ancora piuttosto basso, suggerisce che la separazione tra i cluster non è ottimale, ma non peggiora rispetto ai casi precedenti.

  • Corrected Rand Index pari a 0.9406, valore molto alto, indica che il clustering ottenuto con DIANA è molto simile alla partizione ideale.

  • Variational Information (VI) pari a 0.1541, valore basso, indica una buona corrispondenza con la partizione di riferimento.

Applicare DIANA solo sulle variabili bottom e diagonal migliora il risultato della separazione dei cluster, come evidenziato dal valore più alto di Silhouette (0.6241). Tuttavia, l’indice di Dunn rimane basso, suggerendo che la separazione tra i cluster non è ancora ideale. La concordanza con la partizione ideale (Corrected Rand Index) è molto alta, il che indica che il clustering rispecchia bene una divisione naturale del dataset.

Visualizziamo come il dataset è stato diviso in cluster.

final_diana_clustering_bd <- eclust(banknote_bd_scaled, "diana", stand = FALSE, k = 2, hc_metric = "euclidean")

fviz_cluster(final_diana_clustering_bd, geom = "point", ellipse.type = "norm", palette="jco", ggtheme = theme_minimal())

Vantaggi e Svantaggi del Clustering Gerarchico

TODO #### Vantaggi:

Non richiede la specifica del numero di cluster a priori.

Produce una rappresentazione gerarchica che può essere utile per comprendere le relazioni tra i dati.

Adatto a dataset di piccole e medie dimensioni.

Svantaggi:

TODO

Sensibile alla scelta della misura di distanza e del metodo di linkage.

Può essere computazionalmente costoso per dataset di grandi dimensioni. (quello divisivo)

Una volta effettuata una fusione, non è possibile separare i cluster, rendendo l'algoritmo meno flessibile in caso di errori.

3.2 Clustering Partizionale: K-means e K-medoids (PAM)

3.2.1 K-Means

3.2.1.1 Esecuzione dell’algoritmo con tutte le variabili

In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-means, dopo aver scalato le variabili numeriche del dataset.

library(factoextra)
df <- scale(banknote[, -1])

Successivamente si determina il numero ottimale di cluster da utilizzare nell’algoritmo, impiegando il metodo Elbow.

fviz_nbclust(df, kmeans, method = "wss") +
  geom_vline(xintercept = 4, linetype = 2)+
  labs(subtitle = "Elbow method")

Il metodo Silhouette.

fviz_nbclust(df, kmeans, method = "silhouette")+
  labs(subtitle = "Silhouette method")

Ed infine la GAP analysis:

set.seed(123)
fviz_nbclust(df, kmeans, nstart = 25, method = "gap_stat", nboot = 50)+
  labs(subtitle = "Gap statistic method")

Il primo metodo indica quattro come numero ottimale di cluster, il secondo ne suggerisce due mentre la gap statistic tre. Poiché il dataset è composto da 200 osservazioni, equamente suddivise tra banconote autentiche e contraffatte, si sceglie di impostare il numero di cluster a due e si procede con l’esecuzione dell’algoritmo.

Per prima cosa si osserva una matrice di grafici a dispersione, utile per analizzare le relazioni tra le variabili numeriche del dataset.

set.seed(123)
km.res <- kmeans(df, 2, nstart = 25)
print(km.res)
## K-means clustering with 2 clusters of sizes 108, 92
## 
## Cluster means:
##       Length       Left      Right     Bottom        Top   Diagonal
## 1 -0.1221604  0.5482838  0.6163644  0.6788357  0.5494346 -0.8005001
## 2  0.1434057 -0.6436375 -0.7235582 -0.7968941 -0.6449884  0.9397176
## 
## Clustering vector:
##   [1] 1 2 2 2 2 1 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## Within cluster sum of squares by cluster:
## [1] 396.3549 304.8505
##  (between_SS / total_SS =  41.3 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
aggregate(banknote, by=list(cluster=km.res$cluster), mean)
##   cluster Status Length     Left    Right    Bottom      Top Diagonal
## 1       1     NA 214.85 130.3194 130.2056 10.398148 11.09167 139.5611
## 2       2     NA 214.95 129.8891 129.6641  8.266304 10.13261 141.5663
dd <- cbind(banknote, cluster = km.res$cluster)
head(dd)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       1
km.res$cluster
##   [1] 1 2 2 2 2 1 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
km.res$size
## [1] 108  92
km.res$centers
##       Length       Left      Right     Bottom        Top   Diagonal
## 1 -0.1221604  0.5482838  0.6163644  0.6788357  0.5494346 -0.8005001
## 2  0.1434057 -0.6436375 -0.7235582 -0.7968941 -0.6449884  0.9397176
cl <- km.res$cluster
pairs(df, gap=0, pch=cl, col=c("red", "blue")[cl])

TODO: commentare

Ora si analizza la distribuzione dei due cluster.

fviz_cluster(km.res, 
             data = df,
             palette = c("red", "blue"),
             ellipse.type = "euclid", 
             star.plot = TRUE,
             repel = TRUE,
             ggtheme = theme_minimal()
)

TODO: aggiungere eclust() (che utilizzare gap statistics per calcolare k)

3.2.1.2 Esecuzione dell’algoritmo con solamente le variabili bottom e diagonal

In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-means, dopo aver scalato le variabili numeriche del dataset.

df_ridotto <- df[, c("Bottom", "Diagonal")]

TODO: Successivamente si determina il numero ottimale di cluster da utilizzare nell’algoritmo, impiegando il metodo Elbow.

fviz_nbclust(df_ridotto, kmeans, method = "wss") +
  geom_vline(xintercept = 4, linetype = 2)+
  labs(subtitle = "Elbow method")

TODO: Il metodo Silhouette.

fviz_nbclust(df_ridotto, kmeans, method = "silhouette")+
  labs(subtitle = "Silhouette method")

TODO: Ed infine la GAP analysis:

set.seed(123)
fviz_nbclust(df_ridotto, kmeans, nstart = 25, method = "gap_stat", nboot = 500)+
  labs(subtitle = "Gap statistic method")

TODO: Il primo metodo indica quattro come numero ottimale di cluster, il secondo ne suggerisce due mentre la gap statistic tre. Poiché il dataset è composto da 200 osservazioni, equamente suddivise tra banconote autentiche e contraffatte, si sceglie di impostare il numero di cluster a due e si procede con l’esecuzione dell’algoritmo.

TODO: Per prima cosa si osserva una matrice di grafici a dispersione, utile per analizzare le relazioni tra le variabili numeriche del dataset.

set.seed(123)
km.res <- kmeans(df_ridotto, 2, nstart = 25)
print(km.res)
## K-means clustering with 2 clusters of sizes 99, 101
## 
## Cluster means:
##       Bottom   Diagonal
## 1  0.7829035 -0.9035032
## 2 -0.7674005  0.8856120
## 
## Clustering vector:
##   [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## Within cluster sum of squares by cluster:
## [1] 81.94094 35.86830
##  (between_SS / total_SS =  70.4 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
aggregate(banknote, by=list(cluster=km.res$cluster), mean)
##   cluster Status   Length     Left    Right    Bottom      Top Diagonal
## 1       1     NA 214.8222 130.3000 130.1939 10.548485 11.12727 139.4424
## 2       2     NA 214.9683 129.9465 129.7238  8.308911 10.18317 141.5040
dd <- cbind(banknote, cluster = km.res$cluster)
head(dd)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       2
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       2
km.res$cluster
##   [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
km.res$size
## [1]  99 101
km.res$centers
##       Bottom   Diagonal
## 1  0.7829035 -0.9035032
## 2 -0.7674005  0.8856120
cl <- km.res$cluster
pairs(df_ridotto, gap=0, pch=cl, col=c("red", "blue")[cl])

TODO: commentare

TODO: Ora si analizza la distribuzione dei due cluster.

fviz_cluster(km.res, 
             data = df_ridotto,
             palette = c("red", "blue"),
             ellipse.type = "euclid", 
             star.plot = TRUE,
             repel = TRUE,
             ggtheme = theme_minimal()
)

3.2.1.3 Esecuzione dell’algoritmo con solamente le variabili scelte tramite la PCA

In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-means, dopo aver scalato le variabili numeriche del dataset.

df_pca <- scale(banknote_pca[, -1])

TODO: Successivamente si determina il numero ottimale di cluster da utilizzare nell’algoritmo, impiegando il metodo Elbow.

fviz_nbclust(df_pca, kmeans, method = "wss") +
  geom_vline(xintercept = 4, linetype = 2)+
  labs(subtitle = "Elbow method")

TODO: Il metodo Silhouette.

fviz_nbclust(df_pca, kmeans, method = "silhouette")+
  labs(subtitle = "Silhouette method")

TODO: Ed infine la GAP analysis:

set.seed(123)
fviz_nbclust(df_pca, kmeans, nstart = 25, method = "gap_stat", nboot = 500)+
  labs(subtitle = "Gap statistic method")

TODO: Il primo metodo indica quattro come numero ottimale di cluster, il secondo ne suggerisce due mentre la gap statistic tre. Poiché il dataset è composto da 200 osservazioni, equamente suddivise tra banconote autentiche e contraffatte, si sceglie di impostare il numero di cluster a due e si procede con l’esecuzione dell’algoritmo.

TODO: Per prima cosa si osserva una matrice di grafici a dispersione, utile per analizzare le relazioni tra le variabili numeriche del dataset.

set.seed(123)
km.res <- kmeans(df_pca, 2, nstart = 25)
print(km.res)
## K-means clustering with 2 clusters of sizes 100, 100
## 
## Cluster means:
##         PCA1       PCA2        PCA3
## 1 -0.8574889  0.3668574 -0.02482513
## 2  0.8574889 -0.3668574  0.02482513
## 
## Clustering vector:
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   2   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   1   1   1   1   1   1   1   1   1   2   2   1   1   1   1   1   1   1   1   1 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2   2 
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
##   1   2   2   2   2   2   1   2   2   2   2   2   2   2   2   2   2   2   2   2 
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
## 
## Within cluster sum of squares by cluster:
## [1] 221.1989 201.7035
##  (between_SS / total_SS =  29.2 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
aggregate(banknote, by=list(cluster=km.res$cluster), mean)
##   cluster Status  Length    Left   Right Bottom    Top Diagonal
## 1       1     NA 215.008 129.945 129.713  8.307 10.183  141.466
## 2       2     NA 214.784 130.298 130.200 10.528 11.118  139.501
dd <- cbind(banknote, cluster = km.res$cluster)
head(dd)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       2
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       1
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       1
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       1
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       1
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       1
km.res$cluster
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   2   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   1   1   1   1   1   1   1   1   1   2   2   1   1   1   1   1   1   1   1   1 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2   2 
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
##   1   2   2   2   2   2   1   2   2   2   2   2   2   2   2   2   2   2   2   2 
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2
km.res$size
## [1] 100 100
km.res$centers
##         PCA1       PCA2        PCA3
## 1 -0.8574889  0.3668574 -0.02482513
## 2  0.8574889 -0.3668574  0.02482513
cl <- km.res$cluster
pairs(df_ridotto, gap=0, pch=cl, col=c("red", "blue")[cl])

TODO: commentare

TODO: Ora si analizza la distribuzione dei due cluster.

fviz_cluster(km.res, 
             data = df_pca,
             palette = c("red", "blue"),
             ellipse.type = "euclid", 
             star.plot = TRUE,
             repel = TRUE,
             ggtheme = theme_minimal()
)

3.2.2 K-Medoids

3.2.2.1 Esecuzione dell’algoritmo sul dataset selezionando tutte le variabili

In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-medoids con k impostato a 2.

Come in precedenza, si analizzano prima i grafici a dispersione e successivamente i due cluster.

library(cluster)
head(df, n = 3)
##          Length      Left     Right     Bottom       Top  Diagonal
## [1,] -0.2549435  2.433346  2.829942 -0.2890067 -1.183765 0.4482473
## [2,] -0.7860757 -1.167507 -0.634788 -0.9120152 -1.432847 1.0557460
## [3,] -0.2549435 -1.167507 -0.634788 -0.4966762 -1.308306 1.4896737
pam.res <- pam(df, 2)
print(pam.res)
## Medoids:
##       ID     Length       Left      Right     Bottom        Top  Diagonal
## [1,] 193 -0.5205096  0.4944248  0.6026155  0.9570103  0.5598130 -1.113892
## [2,]  47 -0.2549435 -0.6135300 -0.6347880 -0.7735689 -0.5610584  0.882175
## Clustering vector:
##   [1] 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## Objective function:
##    build     swap 
## 1.842553 1.788867 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       2
## 7 genuine  215.5 129.5 129.7    7.9  9.6    141.6       2
## 8 genuine  214.5 129.6 129.2    7.2 10.7    141.7       2
pam.res$medoids
##          Length       Left      Right     Bottom        Top  Diagonal
## [1,] -0.5205096  0.4944248  0.6026155  0.9570103  0.5598130 -1.113892
## [2,] -0.2549435 -0.6135300 -0.6347880 -0.7735689 -0.5610584  0.882175
head(pam.res$clustering)
## [1] 1 2 2 2 2 2
cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])

Come previsto, i cluster risultano molto simili.

fviz_cluster(pam.res,
             palette = c("#00AFBB", "#FC4E07"),
             ellipse.type = "t",
             repel = TRUE,
             ggtheme = theme_classic()
)

3.2.2.2 Esecuzione dell’algoritmo con solamente le variabili bottom e diagonal

df_ridotto <- df[, c("Bottom", "Diagonal")]
pam.res <- pam(df_ridotto, 2)
print(pam.res)
## Medoids:
##       ID     Bottom   Diagonal
## [1,]  47 -0.7735689  0.8821750
## [2,] 185  0.8877871 -0.8535357
## Clustering vector:
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [149] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [186] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## Objective function:
##     build      swap 
## 0.8831314 0.6414212 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       1
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       1
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       1
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       1
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       1
## 7 genuine  215.5 129.5 129.7    7.9  9.6    141.6       1
## 8 genuine  214.5 129.6 129.2    7.2 10.7    141.7       1
pam.res$medoids
##          Bottom   Diagonal
## [1,] -0.7735689  0.8821750
## [2,]  0.8877871 -0.8535357
head(pam.res$clustering)
## [1] 1 1 1 1 1 1
cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])

Come previsto, i cluster risultano molto simili.

fviz_cluster(pam.res,
             palette = c("#00AFBB", "#FC4E07"),
             ellipse.type = "t",
             repel = TRUE,
             ggtheme = theme_classic()
)

3.2.2.3 Esecuzione dell’algoritmo con solamente le variabili scelte tramite la PCA

df_pca <- scale(banknote_pca[, -1])
pam.res <- pam(df_pca, 2)
print(pam.res)
## Medoids:
##      ID       PCA1       PCA2       PCA3
## 198 198  0.8816755 -0.1289878 0.19323058
## 21   21 -1.1983358  0.0749874 0.05222386
## Clustering vector:
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   2   2   2   2   1   2   2   2   1   1   2   2   2   2   2   2   2   2   2 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2   2   2   2   2   2 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   2   2   2   2   2   1   2   2   2   1   2   2   2   2   2   2   2   2   2   2 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   2   2   2   2   1   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2 
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## Objective function:
##    build     swap 
## 1.370309 1.329691 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       1
## 7 genuine  215.5 129.5 129.7    7.9  9.6    141.6       2
## 8 genuine  214.5 129.6 129.2    7.2 10.7    141.7       2
pam.res$medoids
##           PCA1       PCA2       PCA3
## 198  0.8816755 -0.1289878 0.19323058
## 21  -1.1983358  0.0749874 0.05222386
head(pam.res$clustering)
## 1 2 3 4 5 6 
## 1 2 2 2 2 1
cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])

Come previsto, i cluster risultano molto simili.

fviz_cluster(pam.res,
             palette = c("#00AFBB", "#FC4E07"),
             ellipse.type = "t",
             repel = TRUE,
             ggtheme = theme_classic()
)

3.1.4 Misture di Gaussiane

In questa sezione introduciamo l’applicazione del clustering basato su misture di gaussiane, utilizzando la funzione Mclust del pacchetto mclust. Questo approccio assume che i dati siano generati da una combinazione di distribuzioni gaussiane e, tramite l’algoritmo EM (Expectation-Maximization), stima i parametri di ciascuna componente. Il vantaggio principale di questo metodo è la capacità di modellare cluster con forme ellissoidali, permettendo di gestire anche cluster con sovrapposizioni. Inoltre, Mclust effettua una selezione automatica del numero ottimale di cluster e del modello di covarianza (parsimonioso) in base al criterio BIC, fornendo un’indicazione della bontà dell’adattamento e della complessità del modello.

library(mclust)
data(banknote)
dati <- banknote[, -1]
dati_scaled <- scale(dati)
gmm_result <- Mclust(dati_scaled)
summary(gmm_result)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 4 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -1191.595 200 51 -2653.405 -2666.898
## 
## Clustering table:
##  1  2  3  4 
## 16 99 47 38
plot(gmm_result, what = "BIC")

plot(gmm_result, what = "classification")

plot(gmm_result, what = "uncertainty")

plot(gmm_result, what = "density")

Nonostante sappiamo che il dataset banknote presenta due classi reali (ad esempio, banconote fraudolente e non fraudolente), l’analisi con Mclust ha individuato come miglior modello un numero di cluster pari a 3. Inoltre, il modello selezionato è di tipo VEE (ovvero ellissoidale, con forma variabile, ma con volume ed orientamento uguali), in base al valore più alto di BIC evidenziato nel plot. È interessante notare come il BIC per 2 cluster risulti significativamente inferiore rispetto a quello per un numero maggiore di cluster, probabilmente perché sono state utilizzate tutte le variabili del dataset per costruire i cluster.
Confrontiamo esplicitamente i valori di BIC e ICL per soluzioni con 2 e 3 cluster:

library(mclust)
data(banknote)
dati <- banknote[, -1]
dati_scaled <- scale(dati)
gmm_result2 <- Mclust(dati_scaled, G=2)
summary(gmm_result2)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -1260.823 200 49 -2781.264 -2781.417
## 
## Clustering table:
##   1   2 
## 101  99
gmm_result3 <- Mclust(dati_scaled, G=3)
summary(gmm_result3)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 3 components: 
## 
##  log-likelihood   n df      BIC      ICL
##       -1212.796 200 43 -2653.42 -2653.74
## 
## Clustering table:
##  1  2  3 
## 16 99 85

Dai sommari ottenuti si osserva che entrambi i modelli, per k = 2 e k = 3, forniscono valori che suggeriscono una soluzione migliore per 3 cluster. I valori del BIC e, conseguentemente, anche quelli dell’ICL (che include una penalizzazione per l’incertezza nelle assegnazioni), indicano che la soluzione a 3 cluster è preferibile rispetto a quella a 2 cluster.
Valutiamo l’effetto della scelta delle variabili sul modello. Utilizziamo invece solo le due variabili “Bottom” e “Diagonal”, individuate in precedenza dagli scatterplot come le più discriminanti:

library(mclust)
data(banknote)
dati <- banknote[, c("Bottom", "Diagonal")]
dati_scaled <- scale(dati)
gmm_result <- Mclust(dati_scaled)
summary(gmm_result)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EVE (ellipsoidal, equal volume and orientation) model with 3 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -387.8563 200 13 -844.5908 -849.2887
## 
## Clustering table:
##   1   2   3 
## 100  16  84
plot(gmm_result, what = "BIC")

plot(gmm_result, what = "classification")

plot(gmm_result, what = "uncertainty")

plot(gmm_result, what = "density")

Anche in questo caso si nota che il valore di BIC è notevolmente più alto (migliore) rispetto al caso in cui sono state usate tutte le variabili. L’algoritmo suggerisce ancora un numero ottimale di cluster pari a 3, ma questa volta ha selezionato come modello parsimonioso il modello EVE (ellissoidale, con volume uguale e orientamento variabile). È importante osservare che in uno dei cluster sono presenti solamente 16 osservazioni, il che potrebbe indicare un gruppo di dati che, probabilmente, sono state erroneamente assegnate a un cluster separato.
Verifichiamo ora se i valori di ICL sono in accordo con quelli del BIC confrontando esplicitamente le soluzioni per 2 e 3 cluster:

library(mclust)
data(banknote)
dati <- banknote[, c("Bottom", "Diagonal")]
dati_scaled <- scale(dati)
gmm_result2 <- Mclust(dati_scaled, G=2)
summary(gmm_result2)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVE (ellipsoidal, equal orientation) model with 2 components: 
## 
##  log-likelihood   n df      BIC       ICL
##       -417.7539 200 10 -888.491 -890.1796
## 
## Clustering table:
##   1   2 
##  99 101
gmm_result3 <- Mclust(dati_scaled, G=3)
summary(gmm_result3)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EVE (ellipsoidal, equal volume and orientation) model with 3 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -387.8563 200 13 -844.5908 -849.2887
## 
## Clustering table:
##   1   2   3 
## 100  16  84

In questo caso, i valori di BIC e ICL sono concordi: per un numero di cluster pari a 3 si ottiene un BIC di circa -845 e un ICL di -849, mentre per 2 cluster il BIC risulta attorno a -888 e l’ICL a -890. Questo dimostra che anche l’ICL, che penalizza ulteriormente le assegnazioni incerte, suggerisce che il modello ottimale è quello con 3 cluster.
Adesso, procederemo ad utilizzare solamente le variabili estratte tramite la tecnica PCA.

dati <- banknote_pca[, -1]
dati_scaled <- scale(dati)
gmm_result <- Mclust(dati_scaled)
summary(gmm_result)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -764.5395 200 16 -1613.852 -1616.362
## 
## Clustering table:
##   1   2 
## 102  98
plot(gmm_result, what = "BIC")

plot(gmm_result, what = "classification")

plot(gmm_result, what = "uncertainty")

plot(gmm_result, what = "density")

In questo caso si nota che il valore di BIC è più alto (migliore) rispetto al caso in cui sono state usate tutte le variabili ma è più basso del caso in cui sono state usate solamente le due variabili selezionate. Questa volta l’algoritmo suggerisce un numero ottimale di cluster pari a 2, selezionando come modello parsimonioso il modello EEV (ellissoidale, con volume e orientamento uguale).
Verifichiamo ora se i valori di ICL sono in accordo con quelli del BIC confrontando esplicitamente le soluzioni per 2 e 3 cluster:

gmm_result2 <- Mclust(dati_scaled, G=2)
summary(gmm_result2)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -764.5395 200 16 -1613.852 -1616.362
## 
## Clustering table:
##   1   2 
## 102  98
gmm_result3 <- Mclust(dati_scaled, G=3)
summary(gmm_result3)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 3 components: 
## 
##  log-likelihood   n df       BIC       ICL
##        -764.778 200 19 -1630.224 -1657.822
## 
## Clustering table:
##  1  2  3 
## 64 98 38

Anche in questo caso, i valori di BIC e ICL sono concordi: per un numero di cluster pari a 3 si ottiene un BIC di circa -1630 e un ICL di -1657, mentre per 2 cluster il BIC risulta attorno a -1614 e l’ICL a -1616. Questo dimostra che anche l’ICL suggerisce che il modello ottimale è quello con 2 cluster.
## 3.1.5. Validazione e visualizzazione dei risultati

TODO: inserire commento introduttivo

3.1.5.3. Validazione con la classe reale

In questa sezione analizzeremo la qualità dei diversi metodi di clustering applicati al dataset banknote, fissando il numero di cluster pari al numero di classi reali, ovvero 2. Questo approccio ci permette di utilizzare tecniche di validazione esterne, dato che disponiamo delle etichette reali delle osservazioni (variabile ‘Status’).
Per valutare l’efficacia dei cluster identificati, utilizzeremo diverse metriche, tra cui:
- Adjusted Rand Index (ARI) e Meila Variation Index, che misurano il grado di concordanza tra il clustering ottenuto e la classificazione reale. - Metriche di valutazione proprie del machine learning supervisionato, come accuratezza, specificità, recall, matrice di confusione e AUC, per interpretare il clustering come un problema di classificazione binaria.
Verranno valutate le misure di validazione al variare del numero di variabili utilizzate nel clustering. In particolare, analizzeremo tre configurazioni differenti:
- Dataset completo: utilizzo di tutte le variabili disponibili. - Dataset ridotto: utilizzo esclusivo delle variabili ‘Diagonal’ e ‘Bottom’, identificate come le più discriminative. - Dataset trasformato: utilizzo delle variabili ottenute tramite PCA.

3.1.5.3.1. Dataset completo

Nel seguente codice, applichiamo i diversi algoritmi di clustering al dataset originale con i migliori parametri calcolati in precedenza (fissando il numero di cluster a 2) e aggiungiamo al dataframe le colonne corrispondenti alle assegnazioni di ciascun metodo.

library(cluster) 
data("banknote")
data <- banknote[, -1]
data_scaled <- scale(data)

set.seed(123)
km_model <- kmeans(data_scaled, centers = 2, nstart = 25)
pam_model <- pam(data_scaled, k = 2, nstart = 25)
hc_model <- hclust(dist(data_scaled, method = "euclidean"), method = "ward.D2")
diana_model <- diana(data_scaled)
gmm_model <- Mclust(data_scaled, G = 2)

banknote$kmeans <- km_model$cluster
banknote$kmedoids <- pam_model$clustering
banknote$hierarchical_agglomerative <- cutree(hc_model, k = 2)
banknote$hierarchical_divisive <- cutree(as.hclust(diana_model), k = 2)
banknote$gmm <- gmm_model$classification

head(banknote)
##    Status Length  Left Right Bottom  Top Diagonal kmeans kmedoids
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0      1        1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7      2        2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2      2        2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0      2        2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8      2        2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4      1        2
##   hierarchical_agglomerative hierarchical_divisive gmm
## 1                          1                     1   2
## 2                          2                     2   2
## 3                          2                     2   2
## 4                          2                     2   2
## 5                          2                     2   2
## 6                          2                     2   2

Il codice successivo implementa una funzione per riallineare le assegnazioni dei cluster alle classi reali in modo da rendere i confronti più accurati.

library(caret)
library(pROC)
library(mclust)

bestMapping <- function(true, pred) {

  banknote <- banknote %>%
  mutate(Status_numeric = recode(Status,
                                 "counterfeit" = 1,
                                 "genuine" = 2))


  true <- banknote$Status_numeric
  true <- as.character(true)
  pred <- as.character(pred)
  levels_true <- sort(unique(true))
  
  cm1 <- confusionMatrix(factor(pred, levels = levels_true), factor(true, levels = levels_true))
  acc1 <- cm1$overall["Accuracy"]
  
  pred_swap <- ifelse(pred == levels_true[1], levels_true[2], levels_true[1])
  cm2 <- confusionMatrix(factor(pred_swap, levels = levels_true), factor(true, levels = levels_true))
  acc2 <- cm2$overall["Accuracy"]
  
  if(acc2 > acc1) {
    return(list(mapped = pred_swap, cm = cm2))
  } else {
    return(list(mapped = pred, cm = cm1))
  }
}

Nel blocco successivo, calcoliamo le metriche di valutazione per ciascun algoritmo di clustering. È importante notare che il valore di AUC ha un significato più chiaro per l’algoritmo GMM (Gaussian Mixture Model), poiché questo metodo fornisce direttamente le probabilità a posteriori di appartenenza ai cluster. Per gli altri algoritmi, che operano in modo hard (assegnando ciascun punto a un solo cluster senza una misura di incertezza), il calcolo dell’AUC è stato forzato attraverso:

  • Per K-means e K-medoids, è stato definito uno score basato sulla distanza dei punti dal centroide o dal medoide del rispettivo cluster.
  • Per gli algoritmi gerarchici, non è stato possibile adottare un criterio analogo per costruire uno score continuo, motivo per cui l’AUC non è stato calcolato (indicata come NA nella tabella).
# K-means
ari_kmeans   <- adjustedRandIndex(banknote$Status, banknote$kmeans)
mapping_kmeans <- bestMapping(banknote$Status, banknote$kmeans)
cm_kmeans <- mapping_kmeans$cm
acc_kmeans   <- cm_kmeans$overall["Accuracy"]
sens_kmeans  <- cm_kmeans$byClass["Sensitivity"]
spec_kmeans  <- cm_kmeans$byClass["Specificity"]

distances_km <- t(apply(data_scaled, 1, function(x) {
  apply(km_model$centers, 1, function(cent) sqrt(sum((x - cent)^2)))
}))
score_kmeans <- distances_km[,2] / rowSums(distances_km)
roc_obj_kmeans <- roc(banknote$Status, score_kmeans)
auc_kmeans <- auc(roc_obj_kmeans)

# K-medoids
ari_kmedoids <- adjustedRandIndex(banknote$Status, banknote$kmedoids)
mapping_kmedoids <- bestMapping(banknote$Status, banknote$kmedoids)
cm_kmedoids <- mapping_kmedoids$cm
acc_kmedoids   <- cm_kmedoids$overall["Accuracy"]
sens_kmedoids  <- cm_kmedoids$byClass["Sensitivity"]
spec_kmedoids  <- cm_kmedoids$byClass["Specificity"]

medoids <- as.matrix(pam_model$medoids)
distances_kmed <- t(apply(data_scaled, 1, function(x) {
  apply(medoids, 1, function(med) sqrt(sum((x - med)^2)))
}))
score_kmedoids <- distances_kmed[,2] / rowSums(distances_kmed)
roc_obj_kmedoids <- roc(banknote$Status, score_kmedoids)
auc_kmedoids <- auc(roc_obj_kmedoids)

# Gerarchico Agglomerativo
ari_hierAgg  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_agglomerative)
mapping_hierAgg <- bestMapping(banknote$Status, banknote$hierarchical_agglomerative)
cm_hierAgg <- mapping_hierAgg$cm
acc_hierAgg   <- cm_hierAgg$overall["Accuracy"]
sens_hierAgg  <- cm_hierAgg$byClass["Sensitivity"]
spec_hierAgg  <- cm_hierAgg$byClass["Specificity"]
auc_hierAgg <- NA

# Gerarchico Divisivo
ari_hierDiv  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_divisive)
mapping_hierDiv <- bestMapping(banknote$Status, banknote$hierarchical_divisive)
cm_hierDiv <- mapping_hierDiv$cm
acc_hierDiv   <- cm_hierDiv$overall["Accuracy"]
sens_hierDiv  <- cm_hierDiv$byClass["Sensitivity"]
spec_hierDiv  <- cm_hierDiv$byClass["Specificity"]
auc_hierDiv <- NA

# GMM
ari_gmm <- adjustedRandIndex(banknote$Status, banknote$gmm)
mapping_gmm <- bestMapping(banknote$Status, banknote$gmm)
cm_gmm <- mapping_gmm$cm
acc_gmm   <- cm_gmm$overall["Accuracy"]
sens_gmm  <- cm_gmm$byClass["Sensitivity"]
spec_gmm  <- cm_gmm$byClass["Specificity"]
score_gmm <- gmm_model$z[,2]
roc_obj_gmm <- roc(banknote$Status, score_gmm)
auc_gmm <- auc(roc_obj_gmm)

Infine, raccogliamo i risultati in un dataframe riassuntivo.

results_df <- data.frame(
  Metodo       = c("K-means", "K-medoids", "Gerarchico Agglomerativo", "Gerarchico Divisivo", "GMM"),
  ARI = c(ari_kmeans, ari_kmedoids, ari_hierAgg, ari_hierDiv, ari_gmm),
  Accuracy     = c(acc_kmeans, acc_kmedoids, acc_hierAgg, acc_hierDiv, acc_gmm),
  Sensitivity  = c(sens_kmeans, sens_kmedoids, sens_hierAgg, sens_hierDiv, sens_gmm),
  Specificity  = c(spec_kmeans, spec_kmedoids, spec_hierAgg, spec_hierDiv, spec_gmm),
  AUC          = c(auc_kmeans, auc_kmedoids, auc_hierAgg, auc_hierDiv, auc_gmm)
)
print(results_df)
##                     Metodo       ARI Accuracy Sensitivity Specificity    AUC
## 1                  K-means 0.8456292    0.960        1.00        0.92 0.9973
## 2                K-medoids 0.9406018    0.985        1.00        0.97 0.9991
## 3 Gerarchico Agglomerativo 0.9602001    0.990        1.00        0.98     NA
## 4      Gerarchico Divisivo 0.9406015    0.985        0.99        0.98     NA
## 5                      GMM 0.9799995    0.995        1.00        0.99 0.9999

I risultati mostrano che il metodo GMM ha la migliore performance complessiva, con il più alto valore di ARI (0.98) e un’accuratezza del 99.5%. Il clustering gerarchico agglomerativo ha anch’esso buone prestazioni, mentre K-means e K-medoids mostrano una buona affidabilità ma con AUC leggermente inferiori. Analizzando la sensibilità e la specificità, si nota che la sensibilità – ovvero la capacità di identificare correttamente le banconote contraffatte – è generalmente molto alta per tutti gli algoritmi, indicando un’elevata capacità di individuare i falsi. Tuttavia, la specificità, ossia la capacità di riconoscere correttamente le banconote genuine, varia maggiormente tra i metodi. In particolare, K-means mostra una specificità inferiore (0.92) rispetto a K-medoids (0.97) e ai metodi gerarchici (0.98), suggerendo una maggiore tendenza a classificare erroneamente alcune banconote genuine come contraffatte; questo errore è meno grave di classificare una banconota contraffatta come vera.

3.1.5.3.2. Dataset ridotto

Nel seguente codice, applichiamo i diversi algoritmi di clustering al dataset ridotto (con le sole variabili ‘Diagonal’ e ‘Bottom’) con i migliori parametri calcolati in precedenza (fissando il numero di cluster a 2) e aggiungiamo al dataframe le colonne corrispondenti alle assegnazioni di ciascun metodo.

library(cluster) 
data("banknote")
data <- banknote[, c("Bottom", "Diagonal")]
data_scaled <- scale(data)

set.seed(123)
km_model <- kmeans(data_scaled, centers = 2, nstart = 25)
pam_model <- pam(data_scaled, k = 2, nstart = 25)
#TODO: verificare se ward.D2 è il migliore
hc_model <- hclust(dist(data_scaled, method = "euclidean"), method = "ward.D2")
diana_model <- diana(data_scaled)
gmm_model <- Mclust(data_scaled, G = 2)

banknote$kmeans <- km_model$cluster
banknote$kmedoids <- pam_model$clustering
banknote$hierarchical_agglomerative <- cutree(hc_model, k = 2)
banknote$hierarchical_divisive <- cutree(as.hclust(diana_model), k = 2)
banknote$gmm <- gmm_model$classification

head(banknote)
##    Status Length  Left Right Bottom  Top Diagonal kmeans kmedoids
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0      2        1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7      2        1
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2      2        1
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0      2        1
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8      2        1
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4      2        1
##   hierarchical_agglomerative hierarchical_divisive gmm
## 1                          1                     1   1
## 2                          1                     1   1
## 3                          1                     1   1
## 4                          1                     1   1
## 5                          1                     1   1
## 6                          1                     1   1

Nel blocco successivo, calcoliamo le metriche di valutazione come nel caso precedente.

# K-means
ari_kmeans   <- adjustedRandIndex(banknote$Status, banknote$kmeans)
mapping_kmeans <- bestMapping(banknote$Status, banknote$kmeans)
cm_kmeans <- mapping_kmeans$cm
acc_kmeans   <- cm_kmeans$overall["Accuracy"]
sens_kmeans  <- cm_kmeans$byClass["Sensitivity"]
spec_kmeans  <- cm_kmeans$byClass["Specificity"]

distances_km <- t(apply(data_scaled, 1, function(x) {
  apply(km_model$centers, 1, function(cent) sqrt(sum((x - cent)^2)))
}))
score_kmeans <- distances_km[,2] / rowSums(distances_km)
roc_obj_kmeans <- roc(banknote$Status, score_kmeans)
auc_kmeans <- auc(roc_obj_kmeans)

# K-medoids
ari_kmedoids <- adjustedRandIndex(banknote$Status, banknote$kmedoids)
mapping_kmedoids <- bestMapping(banknote$Status, banknote$kmedoids)
cm_kmedoids <- mapping_kmedoids$cm
acc_kmedoids   <- cm_kmedoids$overall["Accuracy"]
sens_kmedoids  <- cm_kmedoids$byClass["Sensitivity"]
spec_kmedoids  <- cm_kmedoids$byClass["Specificity"]

medoids <- as.matrix(pam_model$medoids)
distances_kmed <- t(apply(data_scaled, 1, function(x) {
  apply(medoids, 1, function(med) sqrt(sum((x - med)^2)))
}))
score_kmedoids <- distances_kmed[,2] / rowSums(distances_kmed)
roc_obj_kmedoids <- roc(banknote$Status, score_kmedoids)
auc_kmedoids <- auc(roc_obj_kmedoids)

# Gerarchico Agglomerativo
ari_hierAgg  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_agglomerative)
mapping_hierAgg <- bestMapping(banknote$Status, banknote$hierarchical_agglomerative)
cm_hierAgg <- mapping_hierAgg$cm
acc_hierAgg   <- cm_hierAgg$overall["Accuracy"]
sens_hierAgg  <- cm_hierAgg$byClass["Sensitivity"]
spec_hierAgg  <- cm_hierAgg$byClass["Specificity"]
auc_hierAgg <- NA

# Gerarchico Divisivo
ari_hierDiv  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_divisive)
mapping_hierDiv <- bestMapping(banknote$Status, banknote$hierarchical_divisive)
cm_hierDiv <- mapping_hierDiv$cm
acc_hierDiv   <- cm_hierDiv$overall["Accuracy"]
sens_hierDiv  <- cm_hierDiv$byClass["Sensitivity"]
spec_hierDiv  <- cm_hierDiv$byClass["Specificity"]
auc_hierDiv <- NA

# GMM
ari_gmm <- adjustedRandIndex(banknote$Status, banknote$gmm)
mapping_gmm <- bestMapping(banknote$Status, banknote$gmm)
cm_gmm <- mapping_gmm$cm
acc_gmm   <- cm_gmm$overall["Accuracy"]
sens_gmm  <- cm_gmm$byClass["Sensitivity"]
spec_gmm  <- cm_gmm$byClass["Specificity"]
score_gmm <- gmm_model$z[,2]
roc_obj_gmm <- roc(banknote$Status, score_gmm)
auc_gmm <- auc(roc_obj_gmm)

Infine, raccogliamo i risultati in un dataframe riassuntivo.

results_df <- data.frame(
  Metodo       = c("K-means", "K-medoids", "Gerarchico Agglomerativo", "Gerarchico Divisivo", "GMM"),
  ARI = c(ari_kmeans, ari_kmedoids, ari_hierAgg, ari_hierDiv, ari_gmm),
  Accuracy     = c(acc_kmeans, acc_kmedoids, acc_hierAgg, acc_hierDiv, acc_gmm),
  Sensitivity  = c(sens_kmeans, sens_kmedoids, sens_hierAgg, sens_hierDiv, sens_gmm),
  Specificity  = c(spec_kmeans, spec_kmedoids, spec_hierAgg, spec_hierDiv, spec_gmm),
  AUC          = c(auc_kmeans, auc_kmedoids, auc_hierAgg, auc_hierDiv, auc_gmm)
)
print(results_df)
##                     Metodo       ARI Accuracy Sensitivity Specificity    AUC
## 1                  K-means 0.9799995    0.995        0.99        1.00 0.9999
## 2                K-medoids 0.9799995    0.995        0.99        1.00 0.9999
## 3 Gerarchico Agglomerativo 0.9602000    0.990        0.99        0.99     NA
## 4      Gerarchico Divisivo 0.9406015    0.985        0.98        0.99     NA
## 5                      GMM 0.9799995    0.995        1.00        0.99 0.9995

Questa configurazione semplifica il modello senza perdere accuratezza, dimostrando che ‘Bottom’ e ‘Diagonal’ sono le variabili chiave per distinguere le classi:
- K-means, K-medoids e GMM ottengono ARI 0.98, accuracy 99.5% e AUC ~1, confermando che queste due variabili sono altamente discriminative.
- GMM raggiunge sensibilità 100%, identificando tutte le banconote contraffatte, mentre K-means e K-medoids garantiscono specificità 100%, evitando falsi positivi.

3.1.5.3.3. Dataset trasformato

Nel seguente codice, applichiamo i diversi algoritmi di clustering al dataset trasformato (con le variabili trasformate mediante tecnica PCA) con i migliori parametri calcolati in precedenza (fissando il numero di cluster a 2) e aggiungiamo al dataframe le colonne corrispondenti alle assegnazioni di ciascun metodo.

data <- banknote_pca[, -1]
data_scaled <- scale(data)

set.seed(123)
km_model <- kmeans(data_scaled, centers = 2, nstart = 25)
pam_model <- pam(data_scaled, k = 2, nstart = 25)
#TODO: verificare se ward.D2 è il migliore
hc_model <- hclust(dist(data_scaled, method = "euclidean"), method = "ward.D2")
diana_model <- diana(data_scaled)
gmm_model <- Mclust(data_scaled, G = 2)

banknote$kmeans <- km_model$cluster
banknote$kmedoids <- pam_model$clustering
banknote$hierarchical_agglomerative <- cutree(hc_model, k = 2)
banknote$hierarchical_divisive <- cutree(as.hclust(diana_model), k = 2)
banknote$gmm <- gmm_model$classification

head(banknote)
##    Status Length  Left Right Bottom  Top Diagonal kmeans kmedoids
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0      2        1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7      1        2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2      1        2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0      1        2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8      1        2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4      1        1
##   hierarchical_agglomerative hierarchical_divisive gmm
## 1                          1                     1   1
## 2                          2                     2   2
## 3                          2                     2   2
## 4                          2                     2   2
## 5                          2                     2   2
## 6                          1                     2   2

Nel blocco successivo, calcoliamo le metriche di valutazione come nei casi precedenti.

# K-means
ari_kmeans   <- adjustedRandIndex(banknote$Status, banknote$kmeans)
mapping_kmeans <- bestMapping(banknote$Status, banknote$kmeans)
cm_kmeans <- mapping_kmeans$cm
acc_kmeans   <- cm_kmeans$overall["Accuracy"]
sens_kmeans  <- cm_kmeans$byClass["Sensitivity"]
spec_kmeans  <- cm_kmeans$byClass["Specificity"]

distances_km <- t(apply(data_scaled, 1, function(x) {
  apply(km_model$centers, 1, function(cent) sqrt(sum((x - cent)^2)))
}))
score_kmeans <- distances_km[,2] / rowSums(distances_km)
roc_obj_kmeans <- roc(banknote$Status, score_kmeans)
auc_kmeans <- auc(roc_obj_kmeans)

# K-medoids
ari_kmedoids <- adjustedRandIndex(banknote$Status, banknote$kmedoids)
mapping_kmedoids <- bestMapping(banknote$Status, banknote$kmedoids)
cm_kmedoids <- mapping_kmedoids$cm
acc_kmedoids   <- cm_kmedoids$overall["Accuracy"]
sens_kmedoids  <- cm_kmedoids$byClass["Sensitivity"]
spec_kmedoids  <- cm_kmedoids$byClass["Specificity"]

medoids <- as.matrix(pam_model$medoids)
distances_kmed <- t(apply(data_scaled, 1, function(x) {
  apply(medoids, 1, function(med) sqrt(sum((x - med)^2)))
}))
score_kmedoids <- distances_kmed[,2] / rowSums(distances_kmed)
roc_obj_kmedoids <- roc(banknote$Status, score_kmedoids)
auc_kmedoids <- auc(roc_obj_kmedoids)

# Gerarchico Agglomerativo
ari_hierAgg  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_agglomerative)
mapping_hierAgg <- bestMapping(banknote$Status, banknote$hierarchical_agglomerative)
cm_hierAgg <- mapping_hierAgg$cm
acc_hierAgg   <- cm_hierAgg$overall["Accuracy"]
sens_hierAgg  <- cm_hierAgg$byClass["Sensitivity"]
spec_hierAgg  <- cm_hierAgg$byClass["Specificity"]
auc_hierAgg <- NA

# Gerarchico Divisivo
ari_hierDiv  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_divisive)
mapping_hierDiv <- bestMapping(banknote$Status, banknote$hierarchical_divisive)
cm_hierDiv <- mapping_hierDiv$cm
acc_hierDiv   <- cm_hierDiv$overall["Accuracy"]
sens_hierDiv  <- cm_hierDiv$byClass["Sensitivity"]
spec_hierDiv  <- cm_hierDiv$byClass["Specificity"]
auc_hierDiv <- NA

# GMM
ari_gmm <- adjustedRandIndex(banknote$Status, banknote$gmm)
mapping_gmm <- bestMapping(banknote$Status, banknote$gmm)
cm_gmm <- mapping_gmm$cm
acc_gmm   <- cm_gmm$overall["Accuracy"]
sens_gmm  <- cm_gmm$byClass["Sensitivity"]
spec_gmm  <- cm_gmm$byClass["Specificity"]
score_gmm <- gmm_model$z[,2]
roc_obj_gmm <- roc(banknote$Status, score_gmm)
auc_gmm <- auc(roc_obj_gmm)

Infine, raccogliamo i risultati in un dataframe riassuntivo.

results_df <- data.frame(
  Metodo       = c("K-means", "K-medoids", "Gerarchico Agglomerativo", "Gerarchico Divisivo", "GMM"),
  ARI = c(ari_kmeans, ari_kmedoids, ari_hierAgg, ari_hierDiv, ari_gmm),
  Accuracy     = c(acc_kmeans, acc_kmedoids, acc_hierAgg, acc_hierDiv, acc_gmm),
  Sensitivity  = c(sens_kmeans, sens_kmedoids, sens_hierAgg, sens_hierDiv, sens_gmm),
  Specificity  = c(spec_kmeans, spec_kmedoids, spec_hierAgg, spec_hierDiv, spec_gmm),
  AUC          = c(auc_kmeans, auc_kmedoids, auc_hierAgg, auc_hierDiv, auc_gmm)
)
print(results_df)
##                     Metodo       ARI Accuracy Sensitivity Specificity    AUC
## 1                  K-means 0.8830121    0.970        0.97        0.97 0.9987
## 2                K-medoids 0.8272389    0.955        1.00        0.91 0.9946
## 3 Gerarchico Agglomerativo 0.5160469    0.860        0.95        0.77     NA
## 4      Gerarchico Divisivo 0.9212040    0.980        0.98        0.98     NA
## 5                      GMM 0.9212042    0.980        0.99        0.97 0.9983

Rispetto ai risultati con il dataset completo o ridotto, l’uso di solo 3 variabili ottenute con PCA ha ridotto la performance globale, in particolare per K-medoids e Gerarchico Agglomerativo. Sebbene l’accuratezza resti elevata, l’ARI e la specificità mostrano una leggera diminuzione, il che suggerisce che la riduzione delle variabili ha compromesso la capacità di questi algoritmi di catturare la struttura dei dati. In particolare, l’algoritmo Gerarchico Agglomerativo mostra una performance molto più bassa in termini di ARI e accuratezza, indicando che l’algoritmo è più sensibile alla perdita di variabili discriminanti.

In conclusione, dal punto di vista della misurazione delle performance con tecniche esterne (avendo a disposizione la vera classe delle osservazioni), l’algoritmo GMM è risultato il migliore nei tre diversi dataset, anche se, tutto sommato, gli altri algoritmi non hanno ottenuto dei cattivi risultati.